load("nelson_data.RData")
attach(kkk)

#Free Speech and Public Order Importance
a <- lm(kspeech ~ vidcdum)
b <- lm(idisord ~ vidcdum)
c <- lm(kspeech ~ vidcdum + idisord)

#Mediation Effect
med.eff <- b$coef[2] * c$coef[3]
#Direct Effect
d.eff <- c$coef[2]
#Total Effect
tot.eff <- med.eff + d.eff

#Mediation Effect Variance Estimate - Sobel
var.1 <-b$coef[2]^2*vcov(c)[3,3] + c$coef[3]^2*vcov(b)[2,2] 

#Total Effect Variance Estimate
var.2 <- vcov(c)[2,2] + var.1 + 2*b$coef[2]*vcov(c)[2,3]

ci.1 <- c(med.eff - qnorm(.975)*sqrt(var.1), med.eff + qnorm(.975)*sqrt(var.1))
ci.2 <- c(tot.eff - qnorm(.975)*sqrt(var.2), tot.eff + qnorm(.975)*sqrt(var.2))

#Effect Without No Interaction Assumption
d <- lm(kspeech ~ vidcdum + idisord + vidcdum:idisord)

#Mediation Effects
delta.0 <- b$coef[2]*d$coef[3] 
delta.1 <- b$coef[2]*(d$coef[3] + d$coef[4])

#Direct Effects
d.eff.0 <- d$coef[2] + d$coef[4]*b$coef[1]
d.eff.1 <- d$coef[2] + d$coef[4]*(b$coef[1] + b$coef[2])

#Total Effect
tot.eff.2 <- delta.0 + d$coef[2] + d$coef[4]*(b$coef[1] + b$coef[2])

#Mediation Effect Variance Estimates
var.d0 <- d$coef[3]^2*vcov(b)[2,2] + b$coef[2]^2*vcov(d)[3,3]
var.d1 <- (d$coef[3] + d$coef[4])^2*vcov(b)[2,2] + b$coef[2]^2*(vcov(d)[3,3] + vcov(d)[4,4] + 2*vcov(d)[4,3])

#Direct Effect Variance Estimates
var.nd.0 <- (vcov(d)[2,2] + b$coef[1]^2*vcov(d)[4,4] + 2*b$coef[1]*vcov(d)[2,4] + d$coef[4]^2*vcov(b)[1,1])
var.nd.1 <- (vcov(d)[2,2] + (b$coef[1] + b$coef[2])^2*vcov(d)[4,4] + 2*(b$coef[1] + b$coef[2])*vcov(d)[2,4] + 
    d$coef[4]^2*(vcov(b)[1,1] + vcov(b)[2,2] + 2*vcov(b)[1,2]))

#Total Effect Variance Estimate
var.tot <- ((d$coef[3] + d$coef[4])^2*vcov(b)[2,2] + b$coef[2]^2*vcov(d)[3,3] + vcov(d)[2,2] + 
    (b$coef[1] + b$coef[2])^2*vcov(d)[4,4] + d$coef[4]^2*vcov(b)[1,1] + 
    2*(d$coef[4]*(d$coef[3] + d$coef[4])*vcov(b)[1,2] + b$coef[2]*vcov(d)[2,3] + 
    (b$coef[1] + b$coef[2])*vcov(d)[2,4] + b$coef[2]*(b$coef[1] + b$coef[2])*vcov(d)[3,4]))


ci.d0 <- c(delta.0 - qnorm(.975)*sqrt(var.d0), delta.0 + qnorm(.975)*sqrt(var.d0))
ci.d1 <- c(delta.1 - qnorm(.975)*sqrt(var.d1), delta.1 + qnorm(.975)*sqrt(var.d1))
ci.nd.0 <- c(d.eff.0 - qnorm(.975)*sqrt(var.nd.0), d.eff.0 + qnorm(.975)*sqrt(var.nd.0))
ci.nd.1 <- c(d.eff.1 - qnorm(.975)*sqrt(var.nd.1), d.eff.1 + qnorm(.975)*sqrt(var.nd.1))
ci.tot <- c(tot.eff.2 - qnorm(.975)*sqrt(var.tot), tot.eff.2 + qnorm(.975)*sqrt(var.tot))

cat("______________________________________ \n")
cat("Result With No-Interaction Assumptions \n")
cat("______________________________________ \n")
cat("Speech Tolerance Average Mediation Effect for Public Order: ", med.eff, "and Sobel 95% CI: ", ci.1, "\n")
cat("Direct Effect of Treatment:", d.eff, "and 95% CI: ", confint(c, "vidcdum"), "\n")
cat("Total Effect: ", tot.eff, "and 95% CI: ", ci.2, "\n\n")

cat("______________________________________ \n")
cat("Result Without No-Interaction Assumptions \n")
cat("______________________________________ \n")
cat("Speech Tolerance Average Mediation Effect For Control: ", delta.0, "and 95% CI: ", ci.d0, "\n")
cat("Speech Tolerance Average Mediation Effect For Treated: ", delta.1, "and 95% CI: ", ci.d1, "\n")
cat("Direct Effect of Control:", d.eff.0, "and 95% CI: ", ci.nd.0, "\n")
cat("Direct Effect of Treated:", d.eff.1, "and 95% CI: ", ci.nd.1, "\n")
cat("Total Effect: ", tot.eff.2, "and 95% CI: ", ci.tot, "\n\n")

